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Abstract. Explicit numerical methods based on Lax-Friedrichs and Leap- Frog finite difference 
approximations are constructed to find the numerical solution of the first-order hyperbolic partial 
differential equation with point-wise delay or advance, i.e., shift in space. The differential 
equation involving point-wise delay and advance models the distribution of the time intervals 
between successive neuronal firings. In this paper, we continue the numerical study which 
was initiated in [T]. We construct higher order numerical approximations and discuss their 
consistency, stability and convergence. The numerical approximations constructed in this paper 
are consistent, stable under CFL condition, and convergent. We also extend our methods to the 
higher space dimensions. Some test examples are included to illustrate our approach. These 
examples verify the theoretical estimates and shows the effect of point-wise delay on the solution. 
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1. Introduction 

Partial differential difference equations or more generally partial functional differential equa- 
tions can be found in several mathematical models of control theory, mathematical biology, 
climate models, mathematical economics, the theory of systems which communicate through 
lossless channels, meteorology, and many other areas, see [21 [3l HI O |6] . Several biological phe- 
nomena can be mathematically modeled by using the time dependent first-order partial differ- 
ential equations of hyperbolic type which contains point-wise delay or shift in space. A detailed 
mathematical analysis and numerical computations of the ordinary delay differential equations 
is discussed by Bellen and Zennaro [7] and references therein. In this paper, we continue the 
study which was initiated in i.e., we consider the following first-order hyperbolic partial 
differential equation having point-wise delay with an initial data uq on domain := (0,X). In 
general it reads 

ut + aux = bu{x — a,t), x G O, t > 
u{x,0) = uo{x), 

u{s,t) = 4>{s,t), V s G [-a,0], t > 0, for a > 
u{X,t) = ipit), for a < 0, 

where a = a(x, t) and b = b{x, t) are sufficiently smooth functions of x and t, a is the value of the 
point-wise delay which is non-zero fixed real number. We also assume that a and b are bounded 
in the entire domain and a is not changing its sign in the entire domain. Let |a(x, t)| < A and 
< B, y{x,t). The unknown function u is defined in the underlying domain and also in 
the intervals [—a, 0] due to the presence of point-wise delay. So our domain is [—a, 0] U [0, X] 
and t > 0. The coefficients are sufficiently smooth functions in these intervals and the unknown 
function u is as smooth as the initial data. Due to the presence of point-wise delay in equation 
([1]), we need a boundary-interval condition in the left side of domain, i.e., in the interval [—a, 0]. 
The equation ([1]) is first-order hyperbolic with difference terms, so it requires one boundary 
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condition according to the direction of characteristics [8j. A similar equation to ([T|) having 
positive shift in right side is obtained by replacing —a with /3, see pj. Numerical schemes for 
equation having positive shift can be obtained in similar fashion by replacing point-wise delay 
with advance. 

Due to the presence of point-wise delay and non-constant coefficients, it is not difficult but to 
some extent impossible to find the analytic solution of such type of partial differential equations 
by using the usual methods to find the exact solution of partial differential equations, see (9j. 
Also we cannot solve such type of partial differential difference equations with classical numer- 
ical methods |10l W\\ [121 [13] . Ordinary differential equations with difference terms (delay) are 
quite well understood by now but there is no comparable theory for partial differential equations 
(i.e. for time and space dependent unknowns). Numerical solution of a general class of delay 
differential equation, including stiff problem, differential-algebraic delay equations and neutral 
problem is discussed by Guglielmi et al. in [T3]. Implicit Runge-Kutta method is applied in 
modified form and possible difficulties are discussed. Asymptotic stabilities properties of im- 
plicit Runge-Kutta method for ordinary delay differential equations is considered by Hairer et 
al. in |15] . Parabolic partial differential equations with time delay is considered in [16]. The 
equation ([1]) is of the form 1-D scalar conservation laws with source term in right side which is 
the function of x and t and this source contains point-wise delay. 

If delay and advance arguments are sufficiently small, the authors used the Taylor series approx- 
imations for the difference arguments and proposed an explicit numerical scheme based on the 
upwind finite difference method which is discussed in [17]. This method has restriction on the 
size of point-wise delay and advance. Equation containing only point-wise delay or advance is 
considered by the authors in [IJ. This scheme tackles small as well as large difference argument 
and has first-order convergence in both space and time. 

In the following, we construct numerical approximations based on the Lax-Friedrichs and the 
Leap-Prog finite difference approximations to increase the order of convergence. The applicabil- 
ity of the second-order finite difference method remains valid to solve such type of hyperbolic 
partial differential difference equations provided the geometry of underlying domain is not com- 
plicated. Also it is easy to extend our ideas into higher space dimensions in finite difference 
approximations. In the construction of numerical schemes, a special type of mesh is generated 
during the discretization so that the difference argument also belong to the discrete set of grid 
points. The numerical methods works for large as well as small values of point-wise delay. We 
construct the numerical schemes to find the approximate solution of problem ([1]) in section 2 
and discuss the consistency, stability and convergence. Also we discuss the extensions of numer- 
ical methods in higher space dimensions in section 3. In section 4, we include some numerical 
test examples to validate the predicted theory. Finally, in section 5, we make some concluding 
remarks illustrating the effect of difference arguments on the solution behavior. 

2. Numerical Approximations 

In this section, we construct numerical schemes based on the finite difference method [HI [T8| 
119] . We construct explicit numerical approximations for the given equation ([T]) based on Lax- 
Friedrichs and Leap-frog finite difference approximations. For space-time approximations based 
on finite differences, the (x, t) plane is descretized by taking mesh width Ax and time step At, 
and defining the grid points (xj,t„) by 

Xj = jAx, j = 0, 1, J — 1, J; tn = nAt, n = 0, 1, 2.... 

Now we look for the approximate numerical solution Uj' that approximate u{xj,tn),yj,n. We 
write the closure of CIax as J^Ax and = {xj = jAx, j = 0, 1, 2, J). 



FINITE DIFFERENCE APPROXIMATIONS FOR HYPERBOLIC PDF WITH POINT- WISE DELAY 



3 



2.1. Lax-Priedrichs approximation: In this numerical approximation, we approximate the 
time derivative by forward difference and space by centered difference and then we replace C/" by 
the mean value between C^j^i ™^ ^j-i stability purpose, see jB]- Numerical approximation 
is given by 

Jjn \TjrL 

_ jjn _ jjn 

(2) ^ + = b-u{x, - a, 

To tackle the point- wise delay in the numerical approximation ([2]), we discretize the domain 
in such a way that [xj — a) is a nodal point, V j = 0,1,..., J, i.e., we choose Ax s.t. a = 
moAxjiTiQ £ N and we take total number of points in x— direction s.t. 

Ax a 

where mantissa of any real number is defined as positive fractional part of that number. 
The term containing point-wise delay (V j = 0, 1, J) can be written as 

u{xj — a, tn) = u{jAx — UIqAx, tn) 
= u{{j - mo)Ax,tn) 

« [/" 

Therefore, the numerical approximation is given by 

rrn i j rti 

jjn+l _ ^.+1+^.-1 Jjn _ 

(-i) 2 I r,n J + 1 izl — h"'TT'^ V ?' — 1 2 7-1 

^"j 2Ax ~"j^J-"»0' VJ-i,Z,...,J 1, 

together with initial and boundary-interval conditions as following 
(4a) U^ = u\xj), Vi = 1,2,..., J-1 

(4b) C/o" = n(s,t„) = 0(s,t„), VsG [-a,0], n = 0,1,2... 

(4c) uy = ^l;itn), n = 0,l,2.... 

Now we prove consistency, stability, and convergence of the numerical approxima- 
tion (3-4). 

Theorem 2.2. Let the function uq be sufficiently smooth then the explicit numerical scheme ^ 
is consistent of order 2 in space and of order 1 in time. Also if CFL condition is satisfied, then 
on the arbitrary finite time interval [0,tf], the numerical scheme (0) is stable in the maximum 
norm where stability constant is of the form C = 1 + 0{At), and hence convergent. 

Proof. Firstly, we study the consistency of the numerical scheme ([2]). The consistency error 
of the numerical scheme ([3]) is the difference between the two sides of the equation when the 
approximation [/" is replaced throughout by the exact solution u{xj,tn) of the differential equa- 
tion ([T|). If u is sufficiently smooth, then the consistency error Tj^ of the difference scheme ([3]) 
is given by 

rn _ ^(x„WO- -(-->^'*"^+"(-'-^'*") n(x,+i,^„)-n(x,,i,y _ 

^ - At ^ ^ 2Ax "j^'j-mo- 

By using the Taylor series approximation for the term u{xj,tn+i) w.r.t. to t and for the terms 
u{xj-i,tn) and n(xj+i,t„) w.r.t. x, we get 

Af At^ At'^ At^ 

I- = K + -Y^tt - ^u,. + 0(At2) + 0{^)\^ + [a{u, + -^u,,,) + 0{Ax')]^ - b]u]_^^ 

At Ax^ Ax^ 



[ut + au^]] - by]_^^ + ^utt - ^u^x + ^a>xxx + 0(Ax^ + At ^Ax^ + At^) 
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As u satisfy the given differential equation ([1]), we have 

[ut + au,]^" - = 0. 

Therefore, the consistency error is given by 



XXX 
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+ 0(Ax^ + At'^Ax'^ + Af 



Therefore, T" — )• while (Ax, At) (0,0), which implies the numerical scheme is consistent of 
order 2 in space and of order 1 in time as long as At^^Ax^ — > 0. 
Now by solving the finite difference scheme ([3]) for t^"^^, we get 



(5) U^+^ = \{U^^, + U^_,) 



a-At 



2Ax 

By applying the triangle inequality, we obtain 



(c/]Vi - uu) + ^tmm-m,, v i = i, 2, j - 1. 



n+li 



< 



1 



a"At 

1 - 

Ax 



\u: 



+ 



1 



a" At 
Ax 



j-mo- 



Now by taking the maximum norm, we get 



\U' 



n+l I 



00, Ax = ma,x\U, 



< — max 
" 2 j 



n+l I 



a" At 
Ax 



1 + 



a^At 
Ax 



|J7JLi| + Atmaxl^J^I {U^ 



Using CFL condition < 1 [where A is the bound of a(x, t), V (x, t)], first two terms in above 
inequality can be combined and we get 



n+l I 



3,Ax < (l + SAt)||[/"||oo,A. 



where |6(x,t)| < 5, V (x,t). The term BAt can be controlled by At from which we can predict 
that effect of the term i?At is of the form 0(At). Using these values, we get the following 
estimate 

||C/"+loo,Ax < (l + 0(At))|IC/"||oo,Ax 



I.e., 



\U'"'~^^\\oo,Ax ^ C||f^"'||oo,Aa:) 



which implies the stability of the numerical approximation, where stability constant C is of the 
form C = l + 0(At). 

Now the error in the numerical approximation ([5]) is given by 

We set = U'^ — u{xj,tn) in ([5]). The approximate solution [/" satisfies the difference equation 
([5]) exactly, while exact solution u{xj,tn) leaves the remainder T"At. Therefore, the error in 
the numerical approximation is given by 



At 



At 



Ax 



and 



0. 



Let ^" = maxj{|e"|,j = 0, 1,..., J}. 
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For |a"^| < 1, 



11/ — max |e • | 



<E'' + \b]\AtE'' + At max IT/] 

< ^" + BAtE'^ + At max |rf | 

j 

= (l + 5At)^" + Atmax|rf|, 

j 

since we are using the given initial value for [/", so E^ = and if we suppose that the consistency 
error is bounded i.e. |r"| < Tmax, then by using induction method in the above inequality 

E ^ nAtTfi-idx ^ tfTfnaxi 

where nAt = t/, which proves that the numerical scheme ([3]) is convergent provided that the 
solution u has bounded derivatives up to second order. □ □ 

2.3. Leap-Frog approximation: In this numerical approximation, we use central difference 
for both the space and time. The numerical approximation for equation ([T]) is given by 

Tjn+l jjn—l jjn jjn 

(6) +a,- -G,., V j - 1, 2, J - 1, 

together with initial and boundary-interval conditions as given in Q. 
We write the numerical approximation as below 

(7) U^^' = U^~' - ^(f/jVi - Uj-i) + 2At6,"C/jL„„, V j = 1, 2, J - 1. 

This scheme uses two time intervals to get a central time difference and spreads its legs to pick 
up the space difference at the intermediate time level. It is an explicit scheme which requires 
a special technique to get it started. The initial condition will usually determine the values of 
but can be obtained by any convenient one-step scheme. In this case, we initialize the 
numerical scheme ([7]) with the Lax-Friedrichs scheme Then this scheme gives U"^, U^,... in 
succession. 

Now we discuss consistency, stability and convergence of this numerical approximation. The 
consistency error of the numerical scheme ([6]) is the difference between the two sides of the 
equation when the approximation [/" is replaced throughout by the exact solution u{xj,tn) of 
the differential equation ([1]). If n is sufficiently smooth, then the consistency error T" of the 
difference scheme ([6]) is given by 



rpn _ U{xj,tn+l) -U{xj,tn-l) ^ u{xj+l,tn) - u{xj-l,tn) _ „ 

^ - 2At ^ ^' 2Ax ^J^'j-n.o 

= [ut + —uttt + 0{At^)\] + [a(n, + -^u^.^) + ©(Ax^)]" - 



At^ Ax^ 



Here we used the Taylor series approximations for the terms u{xj,tn+i) and u{xj-i,tn) w.r.t. t 
and X respectively. 

As u satisfy the given differential equation ([TJ, 

[ut + au^]] - b]u]_^^ = 0. 



6 



PARAMJEET SINGH AND KAPIL K. SHARMA 



Therefore, the consistency error is given by 

■'0 o 

Now — )• while (Ax, At) (0, 0), which imphes that the numerical scheme ^ is consistent 
of order 2 in both space and and time. 

Now the Leap-Frog scheme ([6]) is consistent of order 2 in both space and time, the stability 
condition is the same as in the case of right hand side is zero, see |19|. When the term containing 
point-wise delay is zero, it has been shown |19j that Leap-Frog scheme is stable if the CFL 
number is strictly less than 1, that is < 1, where A is the bound of a. Also we have 

proved in the previous numerical scheme ([3]) that right hand side effects only on the form of 
stability constant. Therefore, the numerical approximation is stable in the L2— norm provided 
CFL condition is satisfied. Now the proposed finite difference approximation is linear, therefore, 
by Lax-Richtmyer Equivalence Theorem it is convergent [19j. 

3. Extension to Higher Spatial Dimensions 

We now consider the extensions of our numerical schemes to the higher spatial dimensions. 
For the sake of simplicity, we consider the problem in two spatial dimensions. The extension to 
three spatial dimensions can be done in a similar fashion [20\ 121]. The natural generalization of 
the one-dimensional model problem ([1]) is the following equation together with initial data and 
boundary-intervals conditions 

(8) ut + + buy = cu{x — a,y — /3,t), 

where a, b, c are functions of x, y and t. a and /3 are the values of the point-wise delay 
in X and y-direction respectively. We consider the rectangular domain in the (x,?/)— plane as 
0<x<X, < y < Y . For numerical approximations, we discretize the domain by taking 
uniform grid points with a spacing Ax in the x— direction and Ay in the y— direction. The grid 
points {xi,yj,tn) are defined as following 

Xi = iAx, i = 0,1, ...,Jx; yj = jAy, j = 0,1, ...,Jy; t„ = nAt, n = 0, 1, 2.... 

Now we write the extension of Lax-Friedrichs scheme ([3]) . The approximate solution is denoted 
by Ul^j. Numerical scheme is given by 

(Tjn —TT^ _i_ Ay-/-" r/" 

2Ay '^^^ '-^'--^ ^''H,j^i-moj-(?o' 

together with appropriate initial data and boundary-interval conditions. Here we take the grid 
points in both the directions (x and y) in such a way that the term containing point-wise delays 
is also belong to discrete set of grid points which can be done as we did in the one dimensional 
case. We take total number of points in both x and y direction such that corresponding delays 
are on rriQ and node and total number of points in both the directions are given by 

X , _,,mantissa(a) , 

J^ = — = kX keN, 

Ax a 

and J, = ^ = ry "^""^7"^^\ r g N. 
^ Ay /3 ' 

Most of the analysis of the numerical approximation in one dimension is easily extended to the 

two dimensional case as discussed by Morton et al. in [8] . Truncation error of this approximation 
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([9]) will remain as in one dimensional case except some additions due to the presence of y variable 
[8]. Usual analysis will give the CFL condition for stability in the following form 

A/\t BAt ^ ^ 
Ax Ax ~ ' 
where A and B are the bounds of a and b respectively. 

Proof of convergence follows in similar way, leading to error in the approximation 

provided that the CFL condition is satisfied and u has bounded derivatives of second order. 
Similarly the extension of the Leap-Frog numerical approximation ([6]) with appropriate initial 
data and boundary-interval conditions is given by 

a" At 6" At 

(10) u-;' = - ^(^r+i, - ut,,) - ^(f^M'+i - + 2At<,f/r-™o.-.o- 

4. Illustration Examples 

Purpose of this section is to include some test examples to validate the predicted theory 
established in the paper and to illustrate the effect of point-wise delay on the solution behavior. 
The maximum absolute errors for the Lax-Friedrichs approximation and square root errors for 
the Leap-Frog approximation are calculated using the double mesh principle [22] as the exact 
solution for the considered examples are not available. We perform numerical computations 
using MATLAB. The maximum absolute error is given by 

i?(Ax, At) = ^^^.^rnax^^^ |[/fi(i,n) - U^'J^^{2j,2n) 
Similarly we find the square root errors for the Leap-Frog approximation. 



Example 1. We consider the differential difference equation ([T]) with the following coefficients 
and initial— boundary conditions : 

a{x,t) = i_,_2x-t+2x^+a;4 ' 6(x, t) = 0.5; 

n(x, 0) = exp[-10(4x - 1)^]; u{s, t) = 0, V s G [-a, 0]. 

We consider = (0, 1), Ax = .001 = At. The numerical solution is plotted in Figure 1 with 
a = 0.02 at the time t = 0.5 by both the numerical approximations. We observe that graphs 
of computed solution are very close to each other but there is a slight difference near the peak 
of the graphs. In Figure 2 and 3, we show the change in solution with the time by both the 
numerical approximations. As time increases, both the graphs shift to the right side with the 
time. The error tables 1 and 2 illustrate that the proposed numerical methods are convergent 
in both space and time direction. Error Tables are plotted by refining the grid points . 



Example 2. We consider the differential difference equation ([T]) with the following variable 
coefficients and initial— boundary conditions : 

u(x, 0) = exp[-10(4x - 1)2]; u{s, t) = 0, V s € [-a, 0]. 

In this example, we consider the original problem ([1]) with variable coefficients. We consider 
Q, = (0, 1), Ax = .001 = At. The computed approximate solution is plotted in Figure 4 with 
a = 0.05 at time t = 0.5 by both the numerical approximations. Both the approximations has 
slightly different at maxima but at other points graphs are very close to each other. To show 
the effect of point-wise delay on the solution behavior, we show the numerical solution with 
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spatial variable x in Figure 5 and 6 by the Lax-Priedrichs and the Leap-Frog approximation 
respectively. By changing the value of point-wise delay, we observe that as the value of a is 
increased, the height of impulse is decreased and width is decreased. As the exact solution of 
this problem is not available, we calculate the maximum and square errors by refining the grid 
points. Analysis of Tables 3 and 4 also verify the convergence of both the approximations in the 
space as well as time. 

Example 3. We consider the 2 — D differential difference equation ([8]) with the following 
coefficients and initial— boundary conditions : 

aix, y, t) = i+2(x+y)t+20r^^+j/^)+x4; ^(^' = ' ^' = ^-^^ 

u{x, 0) = exp[-10(4x + 4y- if]; 

u{si,S2,t) = 0, y si e [-a,0] and Vs2e[-/3,0]. 

We consider the two dimensional problem ([8]) with variable coefficients. We consider = 
(0, 1)X{0, 1), Ax = Ay = .01 and time step At = .001. The approximate numerical solutions 
are plotted with a = 0.5 and (3 = 0.5 at time t = 0.5 by Lax-Friedrichs and Leap-Frog scheme 
in Figure 7 and Figure 8 respectively. 
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Table 1. The maximum absolute errors for Example 1 by using Lax-Priedrichs 
scheme with a = 0.05 



Ati Ax- 



1/100 1/200 1/400 1/800 



Ax/2 0.053623 0.024289 

Ax/4 0.026758 0.011382 

Ax/8 0.014642 0.005163 

Ax/16 0.008525 0.004416 



0.011064 0.004520 

0.006565 0.002265 

0.002344 0.001041 

0.001789 0.000836 



Table 2 . The square root errors for Example 1 by using Leap-Frog scheme with 
a = 0.5 



AtiAx- 



1/100 1/200 1/400 1/800 



Ax/2 0.039607 0.018968 0.009380 0.003678 

Ax/4 0.017185 0.007408 0.003681 0.001538 

Ax/8 0.008517 0.003695 0.001840 0.000769 

Aa:/16 0.003750 0.001347 0.000770 0.000285 

Table 3. The maximum absolute errors for Example 2 by using Lax-Priedrichs 
scheme with a = 0.05 



AtjAx^ 1/100 1/200 1/400 1/800 

Ax/2 0.053308 0.024327 0.011097 0.006541 

Ax/4 0.027900 0.011334 0.005577 0.002174 

Ax/8 0.017055 0.004971 0.002233 0.001044 

Aa;/16 0.013255 0.002274 0.001042 0.000433 



Table 4. The square root errors for Example 2 by using Leap-Frog scheme with 
a = 0.1 



AUAa;^ 1/100 1/200 1/400 1/800 

Aa;/2 0.041308 0.019827 0.009815 0.004897 

Aa;/4 0.020042 0.009843 0.004901 0.002448 

Aa;/8 0.009951 0.004915 0.002450 0.001224 

Aa:/16 0.004969 0.002457 0.001225 0.000612 
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Figure 1. Example 1: The approximate solution by both the schemes for a 
0.8 at t = 0.5 
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Figure 2. Example 1: The approximate solution at different time levels with 
the Lax-Friedrichs scheme for a = 0.08 
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Figure 3. Example 1: The approximate solution at different time levels with 
the Leap-Frog scheme for a = 0.1 
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Figure 4. Example 2: The approximate solution by both the schemes for a = 
.05 at t = 0.5 
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Figure 6. Example 2: The effect of the point- wise delay on solution at t = 0.5 
by Leap-Frog scheme 
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Figure 7. Example 3: The approximate solution for a = 0.5 and /3 = 0.5 at t = 0.5 
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Figure 8. Example 3: The approximate solution for a = 0.5 and /3 = 0.5 at t = 0.5 
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5. Concluding Remarks 

We proposed in this paper a way to process hyperbolic partial differential equation with 
point-wise delay. It mainly relies on the constructions of finite difference approximations of 
order greater than one in space. In this numerical analysis, two explicit numerical methods 
based on the Lax-Friedrichs and the Leap- Frog finite difference methods are constructed to find 
the numerical solution of target problem with point-wise delay. The consistency, stability and 
convergence analysis proves that proposed numerical schemes are consistent, stable with CFL 
condition and convergent in both space and time. These second-order numerical methods in 
space maintains the height and width better than a first-order scheme as discussed by the au- 
thors in [IJ. The effect of point-wise delay on the solution behavior is shown by taking some test 
examples. Error Tables illustrate the fact that methods are convergent in space and time. Also 
we extends our ideas in higher space dimensions and include numerical experiments to show the 
behavior of solution in two space dimensions. 

Finally, we remark that the strategy developed here can be applied to a problem having multiple 
point-wise delay or advance or both. 
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